####################################
# This is code to replicate the analyses and figures from the supplementary materials
# Produces Table 20
# Produces Figure 6
# (Un)Natural Disasters: Electoral Cycles in Disaster Relief
# Code developed by Alicia Cooperman
####################################
# Analysis conducted using
# R version 4.0.2 (2020-06-22) 
# on Platform: x86_64-apple-darwin17.0 (64-bit)
####################################

rm(list=ls(all=TRUE))
dev.off()
library(dplyr) # Version 0.8.5
library(ggplot2) # Version 3.3.0
library(reshape2) # Version 1.4.4
library(xtable) # Version 1.8.4
library(lfe) # Version 2.8.5
load("cooperman_2021_data.Rdata") 
load("cooperman_2021_precip_year.Rdata")


df$ibgecode <- as.character(df$ibgecode)
df$year <- as.character(df$year)

mun_in_dataset <- unique(df$ibgecode)
year_in_dataset <- unique(df$year)
precip_year <- precip_year[precip_year$ibgecode %in% mun_in_dataset,]
precip_year <- precip_year[!(precip_year$year %in% year_in_dataset),]

df2 <- bind_rows(df, precip_year)
df2 <- df2 %>% arrange(ibgecode,year)

elec.may <- felm(drought_bin ~ mayorelecyear*belowave6mo_June + time + time.sq + PT + PET_JantoJune +  lncattle + corn_share + bean_share  |ibgecode|0|state.year, data=df)
summary(elec.may)
ate.obs.mayorelecyear <- coef(elec.may)["mayorelecyear"]
ate.obs.mayorelecyear.below <- coef(elec.may)["mayorelecyear:belowave6mo_June"]
ate.obs.below <- coef(elec.may)["belowave6mo_June"]
c(ate.obs.mayorelecyear, ate.obs.mayorelecyear.below, ate.obs.below)

nsims = 1000


# Permutations where each municipality is independent ---------------------

set.seed(123)
assign_independently <- function(df){
  df_sim<-
    df %>%
    group_by(ibgecode) %>%
    mutate(Z_sim = sample(belowave6mo_June, n(), replace=T)) %>%
    arrange(ibgecode)
  return(df_sim$Z_sim)
}
df2 <- df2 %>% arrange(ibgecode,year)
set.seed(123)
# zindep is a matrix in which the columns correspond to different permutations of the treatment assignment vector generated by the function assign_independently_inch()
zindep <- matrix(NA, nrow=nrow(df2), ncol=nsims)
colnames(zindep) <- paste0("V", 1:nsims)
for (i in 1:nsims){
  zindep[,i] <- assign_independently(df2)
  # print(i) # Uncomment to keep track of progress
}
data.indep <- cbind(df2, zindep)
zsims.cols <- colnames(data.indep)[grep("V",colnames(data.indep),fixed=TRUE)]
# Return the estimated ATE under each permutation.
ate.sims.indep.mayorelecyear.below <- rep(NA, nsims)
ate.sims.indep.mayorelecyear <- rep(NA, nsims)
ate.sims.indep.below <- rep(NA, nsims)
for (i in 1:nsims){
  model <- felm(drought_bin ~ mayorelecyear*get(zsims.cols[i]) + time+ time.sq +PT + PET_JantoJune + lncattle + corn_share + bean_share |ibgecode|0|state.year, data=data.indep)
  ate.sims.indep.mayorelecyear.below[i] <- coef(model)["mayorelecyear:get(zsims.cols[i])"]
  ate.sims.indep.mayorelecyear[i] <- coef(model)["mayorelecyear"]
  ate.sims.indep.below[i] <- coef(model)["get(zsims.cols[i])"]
  # print(i) # Uncomment to keep track of progress
}

# Permutations where each state is independent ---------------------

set.seed(123)
assign_state_year <- function(df,rainmeasure){
  wide_df <- dcast(df, ibgecode + state.cr ~ year , value.var = rainmeasure)
  unique_states <- unique(wide_df$state.cr)
  internal_fun <- function(){
    Z <- rep(NA, nrow(wide_df))
    for(i in 1:length(unique_states)){
      local_df <- filter(wide_df, state.cr==unique_states[i])
      Z[wide_df$state.cr==unique_states[i]] <- local_df[,sample(3:ncol(wide_df), 1)]
    }
    return(Z)
  }
  return(as.vector(replicate(length(unique(df$year)),internal_fun())))
}
df2 <- df2 %>% arrange(year,state.cr, ibgecode)
zstate <- matrix(NA, nrow=nrow(df2), ncol=nsims)
colnames(zstate) <- paste0("V", 1:nsims)
for (i in 1:nsims){
  zstate[,i] <- assign_state_year(df2,"belowave6mo_June")
  # print(i) # Uncomment to keep track of progress
}
df2 <- df2 %>% arrange(year,state.cr, ibgecode)
data.state <- cbind(df2,zstate)

zsims.cols <- colnames(data.state)[grep("V",colnames(data.state),fixed=TRUE)]
ate.sims.state.mayorelecyear.below <- rep(NA, nsims)
ate.sims.state.mayorelecyear <- rep(NA, nsims)
ate.sims.state.below <- rep(NA, nsims)
for (i in 1:nsims){
  model <- felm(drought_bin ~ mayorelecyear*get(zsims.cols[i]) + time + time.sq+PT + PET_JantoJune + lncattle + corn_share + bean_share |ibgecode|0|state.year, data=data.state)
  ate.sims.state.mayorelecyear.below[i] <- coef(model)["mayorelecyear:get(zsims.cols[i])"]
  ate.sims.state.mayorelecyear[i] <- coef(model)["mayorelecyear"]
  ate.sims.state.below[i] <- coef(model)["get(zsims.cols[i])"]
  # print(i) # Uncomment to keep track of progress
}

# Permutations where whole Northeast semi-arid together ---------------------

set.seed(123)
assign_NE_year <- function(df,rainmeasure){
  wide_df <- dcast(df, ibgecode + state.cr ~ year , value.var = rainmeasure)
  internal_fun <- function(){
    Z <- wide_df[,sample(3:ncol(wide_df), 1)]
    return(Z)
  }
  return(as.vector(replicate(length(unique(df$year)),internal_fun())))
}
df2 <- df2 %>% arrange(year,state.cr, ibgecode)
zNE <- matrix(NA, nrow=nrow(df2), ncol=nsims)
colnames(zNE) <- paste0("V", 1:nsims)
for (i in 1:nsims){
  zNE[,i] <- assign_NE_year(df2,"belowave6mo_June")
  # print(i) # Uncomment to keep track of progress
}
df2 <- df2 %>% arrange(year,state.cr, ibgecode)
data.NE <- cbind(df2,zNE)

zsims.cols <- colnames(data.NE)[grep("V",colnames(data.NE),fixed=TRUE)]
ate.sims.NE.mayorelecyear.below <- rep(NA, nsims)
ate.sims.NE.mayorelecyear <- rep(NA, nsims)
ate.sims.NE.below <- rep(NA, nsims)
for (i in 1:nsims){
  model <- felm(drought_bin ~ mayorelecyear*get(zsims.cols[i]) + PET_JantoJune + time+ time.sq + PT+ lncattle + corn_share + bean_share  |ibgecode|0|state.year, data=data.NE)
  ate.sims.NE.mayorelecyear.below[i] <- coef(model)["mayorelecyear:get(zsims.cols[i])"]
  ate.sims.NE.mayorelecyear[i] <- coef(model)["mayorelecyear"]
  ate.sims.NE.below[i] <- coef(model)["get(zsims.cols[i])"]
  # print(i) # Uncomment to keep track of progress
}

ate <- list(ate.sims.indep.mayorelecyear.below=ate.sims.indep.mayorelecyear.below, ate.sims.indep.mayorelecyear=ate.sims.indep.mayorelecyear, ate.sims.indep.below=ate.sims.indep.below,
            ate.sims.state.mayorelecyear.below=ate.sims.state.mayorelecyear.below, ate.sims.state.mayorelecyear=ate.sims.state.mayorelecyear, ate.sims.state.below=ate.sims.state.below,
            ate.sims.NE.mayorelecyear.below=ate.sims.NE.mayorelecyear.below,ate.sims.NE.mayorelecyear=ate.sims.NE.mayorelecyear, ate.sims.NE.below=ate.sims.NE.below)


# Table 20: Randomization Inference Estimates of One-Tail P-Values --------

elec.may <- felm(drought_bin ~ mayorelecyear*belowave6mo_June + time + time.sq + PT + PET_JantoJune +  lncattle + corn_share + bean_share  |ibgecode|0|state.year, data=df)
summary(elec.may)
ate.obs.mayorelecyear <- coef(elec.may)["mayorelecyear"]
ate.obs.mayorelecyear.below <- coef(elec.may)["mayorelecyear:belowave6mo_June"]
ate.obs.below <- coef(elec.may)["belowave6mo_June"]
c(ate.obs.mayorelecyear, ate.obs.mayorelecyear.below, ate.obs.below)

one.tail.indep.mayorelecyear.below <- sum(ate$ate.sims.indep.mayorelecyear.below >= ate.obs.mayorelecyear.below)/length(ate$ate.sims.indep.mayorelecyear.below)
one.tail.indep.mayorelecyear <- sum(ate$ate.sims.indep.mayorelecyear <= ate.obs.mayorelecyear)/length(ate$ate.sims.indep.mayorelecyear)
one.tail.indep.below <- sum(ate$ate.sims.indep.below >= ate.obs.below)/length(ate$ate.sims.indep.below)

one.tail.state.mayorelecyear.below <- sum(ate$ate.sims.state.mayorelecyear.below >= ate.obs.mayorelecyear.below)/length(ate$ate.sims.state.mayorelecyear.below)
one.tail.state.mayorelecyear <- sum(ate$ate.sims.state.mayorelecyear <= ate.obs.mayorelecyear)/length(ate$ate.sims.state.mayorelecyear)
one.tail.state.below <- sum(ate$ate.sims.state.below >= ate.obs.below)/length(ate$ate.sims.state.below)

one.tail.NE.mayorelecyear.below <- sum(ate$ate.sims.NE.mayorelecyear.below >= ate.obs.mayorelecyear.below)/length(ate$ate.sims.NE.mayorelecyear.below)
one.tail.NE.mayorelecyear <- sum(ate$ate.sims.NE.mayorelecyear <= ate.obs.mayorelecyear)/length(ate$ate.sims.NE.mayorelecyear)
one.tail.NE.below <- sum(ate$ate.sims.NE.below >= ate.obs.below)/length(ate$ate.sims.NE.below)

RItable.onetail <- cbind(round(c( ate.obs.below, ate.obs.mayorelecyear, ate.obs.mayorelecyear.below),3), 
                         round(c( one.tail.indep.below, one.tail.indep.mayorelecyear, one.tail.indep.mayorelecyear.below),3),
                         round(c( one.tail.state.below, one.tail.state.mayorelecyear, one.tail.state.mayorelecyear.below),3),
                         round(c( one.tail.NE.below, one.tail.NE.mayorelecyear, one.tail.NE.mayorelecyear.below),3))
rownames(RItable.onetail) <- c( "Below Ave Rainfall", "Mayor Elec Year", "Below Ave Rainfall * Mayor Elec Year")
colnames(RItable.onetail) <- c("Coefficient", "Independent", "State", "Northeast Region")
RItable.onetail[RItable.onetail=="0"]<- "<.001"
RItable.onetail

cat(print.xtable(xtable(RItable.onetail,digits=3), floating=F, include.rownames=T), file="SI_tab20.tex", sep="\n")


# Figure 6: Sampling Distribution of Coefficient Estimates --------

coef <- c(rep("Below Average", nsims*3), rep("Below Average * Mayor Elec. Year", nsims*3))
cluster_lev <- rep(c(rep("Municipality", nsims), rep("State", nsims), rep("NE Region", nsims)), 2)
null_ests <- c(ate$ate.sims.indep.below, ate$ate.sims.state.below, ate$ate.sims.NE.below, 
               ate$ate.sims.indep.mayorelecyear.below, ate$ate.sims.state.mayorelecyear.below, ate$ate.sims.NE.mayorelecyear.below)               
ates <- c(rep(ate.obs.below, nsims*3), rep(ate.obs.mayorelecyear.below, nsims*3))    
df3 <- data.frame(cbind(null_ests, cluster_lev, coef, ates), stringsAsFactors = F)

df3$null_ests <- as.numeric(df3$null_ests)
df3$cluster_lev <- factor(df3$cluster_lev, levels=c("Municipality", "State", "NE Region"))
df3$coef <- as.factor(df3$coef)
df3$ates <- as.numeric(df3$ates)

pdf(file = "SI_fig6.pdf", onefile=FALSE,paper="special",width=8,height=8)
par(mar=c(5.1, 5.1, 4.1, 2.1))
P  <- ggplot(df3, aes(x=null_ests)) + geom_histogram(bins=40)  + facet_grid(cluster_lev~coef, scales="free") 
P <- P + geom_vline(aes(xintercept=ates), colour="red", df3) + theme_bw() + xlab("Coef. Estimates Under Sharp Null Hypothesis") + ylab("Count") 
P + theme(strip.text.y = element_text(size = 12, colour = "black", angle = 270), strip.text.x = element_text(size = 12, colour = "black"),axis.title.y = element_text(size = 12),axis.title.x = element_text(size = 12), axis.text = element_text(size=12))
dev.off()
